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The  need  for  new  materials  with  unique  or  improved  properties 
or  as  substitutes  for  critical  materials  that  are  not  native 
to  the  United  States  is  a  continuing  and  critical  problem 
for  the  Department  of  Defense,  for  U.S.  industry  and  for  the 
U.S.  political  posture  in  general.  The  variety  of  techniques 
for  finding  such  materials  has  been  slowly  growing — ranging 
from  exploiting  nature's  creativity  to  employing  a  nuclear 
explosion  to  produce  transuranic  elements  unknown  in  nature. 

One  technique  of  growing  interest  is  rapid  solidification 
and  the  production  thereby  of  metallic  glasses.  By  this 
means  a  number  of  potentially  important  new  materials  are 
being  fabricated  that  could  have  extensive  future  industrial 
uses.  However,  the  importance  of  developing  new  materials 
is  so  great  that  no  avenue  should  go  unexplored. 

In  this  regard,  there  is  one  technique  that  has  been  almost 
entirely  neglected,  despite  the  fact  that  it  has  been 
extremely  successful  in  producing  valuable  new  materials 
in  the  few  cases  where  it  has  been  applied — this  technique 
involves  the  employment  of  elevated  temperature  and  pressure. 
By  this  means  both  diamond  and  cubic  boron  nitride  have  been 
synthesized  and,  in  fact,  are  the  only  materials  so  formed 
that  are  in  major  commercial  production.  What  is  significant 
about  these  materials  is  that  neither  is  thermodynamically 
stable  under  ambient  atmospheric  conditions — although  the 
transition  probability  to  the  normally  stable  state  is  so 


*Other  materials  that  have  been  produced  by  this  technique 
are  principally  minerals  that  have  been  synthesized  for  geo¬ 
physical  research  objectives,  some  simple  inorganic  compounds 
whose  properties  were  not  measured,  and  a  few  organic 
compounds . 
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small  as  to  be  entirely  negligible.  Thus,  with  the  applica¬ 
tion  of  elevated  temperature  and  pressure  it  is  possible  to 
create  important  industrial  materials  that  cannot  be  created 
(in  significant  quantity)  under  normal  atmospheric  conditions 
and  yet  are  metastable  under  these  latter  conditions. 

What  seems  to  be  completely  lacking  is  a  systematic  effort 
to  exploit  the  temperature/pressure  technique — both  experi¬ 
mentally  and  theoretically — for  the  production  of  new  meta¬ 
stable  materials  for  defense  and  industrial  applications. 

The  intention  of  the  work  reported  here  is  to  assist  in  devel¬ 
oping  a  methodology  for  creating  such  materials.  A  wide  range 
of  materials  are  known  to  undergo  transformations  to  new 
phases  with  new  properties  upon  application  of  pressure. 
However,  these  new  phases  change  back  to  the  original  state 
when  pressure  returns  to  ambient. 

Therefore,  the  basic  problem  that  must  be  solved  in  order 
to  fabricate  useful  materials  at  elevated  pressure  is  to 
develop  techniques  to  preserve  the  material  after  pressure 
is  released.  The  solution  of  this  problem  is  related  inti¬ 
mately  to  the  issue  of  how  to  form  suitable  new  metastable 
materials  in  the  first  place.  The  reason  is  that  if  a  new 
state  is  formed  through  a  succession  of  equilibrium  states 
by  applying  pressure,  then  the  process  is  reversible  (except 
possibly  for  minor  time  effects) .  Thus,  either  the  forma¬ 
tion  or  return  to  ambient  conditions,  or  both,  must  involve 
nonequilibrium  states. 

In  the  rapid  solidification  process,  sudden  cooling  inhibits 
the  orderly  transition  through  a  succession  of  equilibrium 
states.  Thus  the  final  state  (e.g.  a  metallic  glass)  is  in 


fact  metastable  with  a  negligible  probability  of  transition 
to  the  stable  crystalline  state.  Accordingly  with  rapid 
solidification,  it  is  the  return  to  ambient  temperatures 
that  involves  nonequilibrium  states. 

The  direct  analog  of  rapid  solidification  in  the  case  of  ele¬ 
vated  pressure  is  to  form  a  new  material  at  an  elevated  pres¬ 
sure  and  then  by  sudden  release  of  this  pressure  to  attempt 
to  maintain  the  new  material  under  ambient  conditions  albeit 
in  a  metastable  state.  Since  the  rate  of  some  mechanisms 
that  induce  transitions  decreases  with  decreasing  pressure, 
it  is  possible  that  the  pressure  analog  of  rapid  solidifica¬ 
tion  will  result  in  new  materials.  Even  though  this  approach 
looks  promising,  it  has  a  disadvantage.  Sudden  changes  of 
pressure  can  be  deleterious  to  high-pressure  equipment.  Thus, 
unless  special  precautions  are  taken  in  the  original  equipment 
design,  this  approach  does  not  appear  attractive.  Since  most 
high-pressure  equipment  was  not  designed  with  this  type  of 
operation  in  mind,  this  approach  should  probably  not  receive 
primary  emphasis — though  it  should  be  explored  to  an  extent 
consistent  with  safe  operation  of  available  high-pressure 
presses,  or  equipment  specially  designed  for  safe  rapid 
decompression  should  be  developed. 

The  primary  experimental  approach*  should  be  analogous  to 
that  employed  in  diamond  synthesis.  This  approach  is  based 
on  the  observation  that  since  the  new  material  to  be  synthe¬ 
sized  must  be  metastable  under  ambient  atmospheric  conditions, 
the  initial  phase  of  the  material  will,  in  general,  be  meta¬ 
stable  in  the  high-pressure/temperature  domain  where  the  new 
material  is  stable.  For  example,  graphite  is  the  stable  form 

* 

Other  approaches  have  been  successfully  employed  in  specific 
cases  but  these  do  not  have  general  applicability. 


of  carbon  for  pressures  less  than  about  20  kbars,  yet  at 
1  bar  the  rate  of  transition  of  diamond  to  graphite  is  so 
slow  as  to  be  entirely  negligible.  In  fact,  only  by  heating 
diamonds  to  about  800°C  at  one  atmosphere  pressure  does  the 
transition  become  observable.  Reciprocally,  one  can  antici¬ 
pate  (as  has  also  been  verified  experimentally)  that  graphite 
will  not  make  a  transition  to  diamond  simply  by  applying 
pressure  in  excess  of  20  kbars,  that  is,  simply  by  entering 
the  diamond-stable  region  of  the  carbon-phase  diagram.  What 
is  required  to  make  the  transition  is  some  method  for  destroy¬ 
ing  the  graphite  structure  and  allowing  the  carbon  to  reform 
as  diamond,  such  as  dissolution  and  reprecipitation  or  some 
form  of  excitation  such  as  heating — perhaps  even  melting  and 
resolidification.  Both  of  the  latter  techniques  have  been 
used  in  synthesizing  diamond  from  graphite. 

More  generally  however,  both  techniques  will  be  necessary 
to  explore  fully  the  phase  diagram  of  materials  at  elevated 
pressure.  If  the  new  phase  is  contiguous  to  the  liquidus 
domain,  heating  perhaps  to  melting  followed  by  cooling  and 
resolidification  can  be  employed.  If  the  new  phase  is  not 
contiguous  to  the  liquidus,  then  dissolution  and  reprecipi¬ 
tation  will  probably  be  necessary. 

In  summary,  the  basic  experimental  approach  should  involve 
the  application  of  increasing  pressure  and,  at  each  pressure, 
either  the  melting  and  resolidification  of  the  specific  mate¬ 
rial  being  studied  or  finding  a  solvent  for  the  material  and 
dissolving  and  reprecipitating  the  material. 

Initially,  it  is  the  choice  of  materials  for  experimental 
investigation  that  presents  the  key  uncertainty  in  exploring 
for  new  metastable  phases.  One  must  be  able  to  define  the 
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region  of  temperature  and  pressure  where  a  particular  phase 
is  stable.  Experimentally  if  one  enters  this  region  through 
a  succession  of  equilibrum  states,  the  entering  phase  will 
in  general  transform  to  the  particular  stable  phase.  On  the 
other  hand,  if  the  entering  phase  is  metastable  relative  to 
the  stable  phase  no  transition  will  take  place.  Experimen¬ 
tally,  therefore  one  would  have  no  indication  that  one  had 
entered  a  region  in  which  another  phase  was  stable  and  there¬ 
fore  that  it  was  worthwhile  attempting  to  convert  to  this 
new  phase.  Only  with  a  priori  knowledge  of  the  phase  diagram 
could  one  be  aware  of  the  new  phase.  To  obtain  this  knowl¬ 
edge,  the  new  phase  must  exist  because  nature  produced  it, 
as  was  the  case  with  diamond,  or  it  must  be  inferred  by  other 
means,  such  as  in  the  case  of  cubic  boron  nitride,  or  more 
generally  one  must  be  able  to  predict  theoretically  the  phase 
diagrams  of  materials.  The  research  described  here  was  an 
attempt  to  develop  a  theoretical  approach  for  guiding  material 
choices  (as  well  as  exploiting  such  past  experimental  data 
as  may  be  available).  Since  the  approach  was  primarily  to 
provide  guidance  for  experiments  rather  than  to  improve  pre¬ 
dictive  accuracy,  emphasis  was  placed  on  developing  simple, 
rapid  and  inexpensive  computational  techniques.  The  present 
report  describes  a  partial  step  in  this  direction. 

Since  the  primary  objective  is  an  approximate  prediction  of 
phase  boundaries,  the  underlying  emphasis  was  placed  on  com¬ 
paring  crystal  structures  for  stability  rather  than  obtaining 
an  accurate  description  of  a  given  structure.  Accordingly, 
those  aspects  of  a  crystal  that  are  insensitive  to  structural 
changes  may  be  neglected  or  only  crudely  estimated. 

The  results  so  far  obtained  apply  to  zero  temperature.  It 
was  not  possible  within  the  scope  of  the  present  effort  to 
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include  effects  of  temperature.  With  this  limitation,  the 
stable  crystal  structure  is  that  one  with  minimum  total 
energy  at  a  given  pressure  or  density.  The  total  energy1 
of  a  crystal  is  generally  expressed  as  a  sum  of  the  ion-ion 
interaction  energy,  the  electron  band  energy,  the  Coulomb 
interaction  energy  of  the  band  electrons,  and  the  band  elec¬ 
tron  exchange/correlation  energy.  The  zero-point  vibrational 
energy  of  the  crystal  and  the  interaction  of  the  crystal  zero 
point  vibration  with  the  electrons  is  usually  ignored,  as  is 
done  here. 

As  a  first  step  in  an  approximate  theory,  only  the  ion-ion 
interaction  energy  and  the  band  energy  have  been  treated 
since  intuitively  the  remaining  terms  which  depend  only  on 
the  electron  density  distribution  appear  less  sensitive  in 
general  to  crystal  structure.  Furthermore,  these  other 
terms  could  add  substantially  to  the  computational  expense. 
Only  if  this  intuitive  judgment  does  not  lead  to  an  approxi¬ 
mate  theory  sufficient  for  guiding  experiment  will  it  be 
necessary  to  develop  approximations  for  the  other  terms. 

The  ion-ion  interaction  energy  is  well-understood  and  is  sum¬ 
marized  in  Section  II.  In  Section  III,  the  band  energy  is 
treated  using  the  so-called  "muffin-tin"  potential.  This 
potential  consists  of  a  spherically  symmetric  potential 
inside  spheres  centered  on  the  ion  sites  and  a  constant 
potential  outside  the  spheres.  Using  the  fact  that  the 
potential  outside  the  spheres  is  a  constant  and  can  be  set 
equal  to  zero,  a  secular  determinant  for  calculating  the 
energy  eigenvalues  is  derived  in  terms  of  the  electron  wave 
function  and  its  derivative  on  the  surface  of  the  spheres. 

The  result  that  we  have  obtained  has  turned  out  not  to  be  new 

2 

but  to  be  closely  related  to  that  of  Korringa  and  equivalent 


to  that  of  Kohn  and  Rostoker,  and  has  been  referred  to  for 
sometime  as  the  KKR  method.  To  carry  out  the  KKR  method  it 
has  been  customary  to  start  with  a  potential,  like  the  Hartree 
Fock,  solve  the  band  problem,  recalculate  the  potential, 
spherically  average  it  inside  the  spheres,  and  calculate  a 
new  average  outside  the  spheres,  re-solve  the  band  problem, 
and  iterate  until  self  consistency  is  achieved.  In  contrast 
we  propose  to  use  a  simple  model  potential  which  if  suffi¬ 
ciently  accurate  to  yield  reasonable  total  band  energies  will 
avoid  the  lengthy  and  costly  iteration.  The  model  potential 
which  we  have  chosen  is  that  of  the  finite  sized  Thomas-Fermi 
(TF)  atom  inside  the  spheres  and  constant  outside.  Using  the 
TF  potential  is  particularly  convenient  since  codes  are  avail¬ 
able  for  calculating  the  surface  values  of  the  wave  function 
and  its  derivative  as  functions  of  atomic  number,  radius  of 
sphere,  energy  and  angular  momentum  quantum  number.  In  addi¬ 
tion,  comparison  of  energy  eigenvalues  using  the  TF  atom  and 
those  derived  for  Hartree-Fock  indicates  comparable  agreement 
with  experiment.* 

In  the  KKR  method  the  most  difficult  quantities  to  calculate 
are  the  so-called  structure  constants.  The  traditional  method 
to  calculate  these  is  based  on  a  technique  due  to  Ewald^  and 
is  the  most  time  consuming  part  of  the  band  energy  calcula¬ 
tion.  Consequently,  we  have  investigated  an  alternate  tech¬ 
nique  that  might  alleviate  this  difficulty.  The  details  of 
the  technique  are  described  in  Section  III.  A  suggestion  for 
possible  further  improvements  is  considered  in  Section  V. 


A  computer  program  that  reflects  the  current  theoretical  model 
has  been  developed.  The  program  for  calculating  the  ion- ion 
interaction  energy  is  discussed  in  Section  II  and  for  the 
band  energy  in  Section  IV. 


In  Section  VI,  we  summarize  what  remains  to  be  done  to 
finalize  the  initial  theoretical  model.  We  also  discuss 
what  would  be  needed  to  validate  the  proposed  approximate 
model  or  to  indicate  that  a  more  complex  model  is  required. 


mm 


II 


ION-ION  INTERACTION 


The  ion-ion  energy  per  unit  cell,  Ecc,  is  defined  as  the 
interaction  energy  of  all  the  ions,  treated  as  point  charge 
in  the  presence  of  a  uniform  background  of  negative  charge 
density  pQ, 


(1) 


The  total  charge  density  is  then 


P(r)  =  e  za  S(z  -  Rn  -  aa)  +  Pq  (2) 

N,o 


and  the  interaction  energy  is  (in  CGS  units) 
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where  the  self  energy  of  the  individual  ions  is  understood 
to  be  excluded  and  X  is  introduced  to  control  the  Coulomb 
singularity.  We  easily  find 


where  the  prime  on  the  sum  indicates  the  N  =  0,  a  -  6  terms 
are  to  be  omitted. 


It  is  convenient  to  remove  the  prime  by  considering  the 
combination 
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an  identity  due  to  Ewald, 
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where  n  is  an  arbitrary  parameter  greater  than  zero.  The 
sensitive  terms  here  are  the  G  =  0  term  in  the  momentum  space 
sum  and  the  N  =  0,  a  =  6  terms  in  the  coordinate  space  sum. 
However,  the  limits  of  these  terms  are  trivial  leading  to  the 
result6 


where  again  the  primes  mean  the  G  =  0  term  and  the  N  =  0, 
a  *  0  terms  are  excluded. 

We  have  developed  a  computer  code  that  calculates  E^.  It 
assumes  all  quantities  are  in  CGS  and  requires  an  input  value 
for  n.  There  is  a  library  routine  that  calculates  the  comple 
mentary  error  function,  erfc(x).  Our  experience  has  been 
that  values  of  n  on  the  order  of  1016  cm-2  are  reasonable 
and  that  the  calculation  is  rapid,  independent  of  the  exact 
value  of  n. 
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III.  BAND  ENERGY 

The  present  method  of  calculating  the  band  energy  is  based  on 
the  "muffin-tin"  approximation  in  which  each  atom  of  the  crys¬ 
tal  is  contained  within  a  sphere  and  the  potential  is  assumed 
constant  outside  the  spheres.  The  underlying  concept  for  the 
method  is  to  determine  the  wave  functions  for  the  valence  and 
conduction  electrons  by  expanding  these  wave  functions  as  a 
superposition  of  atomic  states  within  the  spheres  and  then  to 
match  the  resultant  to  a  superposition  of  eigenstates  external 
to  the  spheres  where  the  potential  is  constant.  This  matching 
is  done  with  the  aid  of  Green's  theorem  relating  the  solution 
external  to  the  spheres  to  the  boundary  values  on  the  spheres. 
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This  method,  while  conceptually  different,  leads  to  identical 

2 

results  to  those  of  Korrmga,  which  is  based  on  scattering 

theory.  An  equivalent  Green's  function  approach  was  developed 

3  5 

by  Kohn  and  Rostoker  and  Ham  and  Segall. 

To  describe  the  method  explicitly,  we  assume  a  general  crys¬ 
tal  with  direct  lattice  vectors  {RN}»  reciprocal  lattice 
vectors  {G }  and  atomic  locations  within  the  unit  cell  [a  }. 
Around  each  atomic  location,  a  sphere  of  radius  Pa  is  circum¬ 
scribed  in  such  a  manner  that  the  spheres  do  not  overlap. 
Within  each  sphere  the  electronic  wave  function  satisfies 


Ha  ib  (p)  =  E  ib  (p)  , 
a  a 


(8) 


where  Ha  is  the  Hamiltonian  of  the  ath  atom  and  E  is 
electron  energy.  Throughout  this  Section,  we  employ 
units  where  distances  are  in  units  of  Bohr  radii  and 
are  in  Rybergs.  The  quantity  kQ  is  defined  as 


the 

natural 

energies 
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k0  =  yrrr  . 


(9) 


The  wave  function,  va  (p  )  ,  can  be  expressed  as  an  expansion  in 
spherical  harmonics, 


a 


i 


,\ 


V1  >  2 

4  ,m 


o 

i  ,  E>0 
1,  ECO 


aL  <  <p>  \m<p>' 


(10) 


where 


h“  R“  (P  )  =  E  Rj  (P  ) 


(ID 


and  for  positive  energies  the  factor  i*  is  introduced  for 

later  convenience.  The  wave  function  in  Eq.  (10)  can  be  con- 

-► 

sidered  to  be  associated  with  the  atom  in  the  R„  =  0  cell; 

o 

wave  functions  in  other  cells  have  an  additional  factor  elk’RN, 
where  k  is  the  electronic  wave  number. 

Outside  the  atomic  spheres,  where  the  potential  is  assumed  to 
be  constant,  the  electron  wave  function,  4»0(r),  can  be 
expressed  in  terms  of  interior  wave  functions  by  means  of 
Green's  theorem.  Thus  we  have 


*  (r)  = 
o 


i  V  f  a  [  ik*RM  -►  3G(r,r 

=  4-  V  Ids  e  N  *  (p  ) - — 

i*  *-•  I  N  I  °  «  3p„ 

N,u/  a 


NaJ 


V  , 


♦  ♦  3\Ji  (o  )  . 

-  elk-R«  G(?'fN«) 
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This  can  be  written  as 
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*0(r> 
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(13) 
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where 


a  ■> 
F  (r) 


ik-RK 


G(r'rNa)' 


-►  ±  -*■ 
r  =  R.T  +  3  +  p  . 

Na  a  a 


Green's  function  satisfies 


(V2  +  E)  G  (r  ,rM  )  =  -4ir6(r-rMJ 


G (r , r  ' )  * 


cosk  r-r' 

_ o' 


,  E  >  0, 


-kQ [ r-r ' 


,  E  <  0. 


Use  has  also  been  made  of  the  equation  satisfied  by  ij*0(r), 

(V2  +  E)  ^Q(r)  =  0,  (18 

since  the  outside  constant  potential  has  been  chosen  as  the 
reference  zero  for  energies. 

The  complete  wave  function  ♦  £(?)  is  thus  determined  by 
V>(r)  =  (p  )  inside  the  ath  sphere,  N  =  0, 


=  ^{r)  outside  the  spheres. 


.  ,  a 

However,  <|£(r)  involves  the  unknown  constants  A^  and  the 

energy  E,  both  of  which  are  functions  of  the  electronic  wave 

number  k.  To  determine  these  unknowns,  one  uses  Eq.  (13) 

and  lets  r  >  a^  +  p^.  Thus  the  requirement 


*o(r>  + 


implies 


CP  3 ) 


*  £s> 
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From  Eq,  (10) ,  it  is  clear  that  Eq.  (21)  constitutes  a  set  of 
linear  equations  in  the  unknowns  A®m.  For  this  set  to  have  a 
nontrivial  solution,  the  determinant  of  the  coefficients  of 
the  Afro's  must  vanish.  Setting  the  determinant  equal  to  zero 
results  in  a  transcendental  equation  for  E  as  a  function  of 
k  (the  so-called  secular  equation).  The  solution  for  E  vs  k 
defines  the  electronic  band  structure  and,  moreover,  deter¬ 
mines  the  values  of  the  A^'s  when  the  requirement  of  wave 
function  normalization  is  added.  Once  E  vs  k  and  the  values 
of  A“m  are  determined,  then  Eq.  (19)  defines  the  overall 
system  wave  function  for  each  k-value.  The  E  vs  k  results 
lead  to  the  total  band  energy  by  summing  to  the  Fermi  energy. 
The  wave  functions  allow  a  calculation  of  the  electronic 
charge  density  which  determines  the  Coulomb  and  correlation 
and  exchange  contributions  to  the  total  energy.  The  most  cri¬ 
tical  element  in  carrying  out  the  above  procedure  is  the 
evaluation  of  Fa(p^  +  a^)  defined  by  Eqs.  (14)  and  (17).  The 
reason  is  that  these  quantities  are  defined  by  infinite 
series  that  are  rather  slowly  convergent  and  may  in  fact 
require  several  thousand  terms  to  achieve  adequate  accuracy. 


Thus  computational  efficiency  and  the  choice  of  form  for 
Fa  (P^+ag )  become  crucial.  Most  of  the  effort  in  exploiting 
the  present  method  has  focused  on  just  these  points. 

The  explicit  form  for  F°(r)  will  now  be  derived.  We  will 
consider  the  positive  energy  case  but  the  corresponding 
equations  for  negative  energy  are  either  obvious  or  will  be 
treated  explicitly.  From  Eqs.  (14)  and  (17)  ,  we  have 


F0(?) 


cosk  r-r„ 
_ o'  Not 


It  is  clear  that  within  a  given  unit  cell  N,  the  only  singu¬ 
larity  of  Fa(r)  results  from  the  single  term 


cosk  | r-r  | 
0 1  Not 


which  satisfies  the  inhomogeneous  Eq.  (16)  .  The  remaining 
infinite  sum  satisfies  the  homogeneous  equation 


-  e^v^NJir.y.!  L  (2J) 

i«.j  / 

Since  the  above  two  properties  apply  for  all  N,  we  can  choose, 
for  convenience,  N  =  0,  for  which 

=  aa  +  ?*'•  <24) 


We  are  interested  in  Fa(r)  when  [cf.  Eq.  (21)] 


and  this  quantity  can  be  considered  to  be  a  function  of  k, 


£  = 


PD  -  P 


B  ct 


(26; 


If  we  define  the  functions  Ha(R)  and  Da  (R)  as 


cos  k-R 


F„  (ao  +  Pfl)  s  Ha(5)  =  6ag  +  Da  (R), 


a  '  g  8 


(27; 


it  is  clear  that  for  sufficiently  small  values  of  R,  Da  (R) 
satisfies  [cf.  Eq.  (23)] 


(?!+E)  Da(S)  =  0. 


(28) 


This  will  be  true  so  long  as  R  is  not  so  large  as  to  encounter 
any  singularities  except  the  possible  one  (for  a  »  e)  at 
R  *  0.  It  is  sufficient  to  assume  that 


R  <  Rn,  N  ^  0,  a  *  8, 


R  <  |%  -  4  +  Rn|  ,  a  ft  0. 


(29) 


Since  Da(R)  satisfies  Eq.  (28),  it  must  have  the  general  form 


°a<R>  “  4*  X^L6 


LM  FAC(L,M)  < 


jL(kQR) ,  E  >  0 
ir  (k_R)  ,  E  <  0 


r*LM<R>' 
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where  ylm<  R)  are  spherical  harmonics,  DLM  are  the  structure 


constants  and  the  factor 


FAC (L ,M)  = 


1 ,  M  >  0 


( -1 )  ,  M  <  0 


2L+1  (L-  M  ) 1 
4 it  (L+  I’m  I)  ! 


has  been  introduced  for  later  convenience.  For  positive 
energies,  jL(x)  are  the  spherical  Bessel  functions  while  for 
negative  energies,  IL(x)  are  the  modified  spherical  Bessel 
functions,  defined  as 

IL(x)  =  (-i) L  jL(ix).  (32; 

Equation  (30)  combined  with  Eq.  (27)  gives  the  most  general 
form  for  Ha(R)  for  sufficiently  small  R. 

An  alternate  form  for  Ha(R)  can  be  derived  by  a  straight¬ 
forward  Fourier  series  transformation  of  Eq.  (22)  multiplied 
by  e-1^*r  since  the  combination  is  then  periodic  in  the  lat¬ 
tice  vectors  RN.  We  easily  find 


H.(J,  .iLSe.i(^S)-.(VV 


(£+3)  2  -  E 


which  is  correct  for  both  positive  and  negative  energies.  By 
making  use  of  the  well-known  relationship, 

ei  (k+G)  *R  „  4ir  £  iL  jL(|ic+G|R)  Y*M(k+G)  Ylm(R),  (34) 
L,M 

we  can  derive  an  explicit  expression  for  the  structure  con- 
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stants,  Dlm,  contained  in 


otB 

dlm  ^ 
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Once  we  have  introduced  the  structure  constants,  we  can 
derive  a  convenient  expression  for  HQ(R)  that  separates  the 
angular  dependences  of  Pg  and  p  This  derivation  employs 
the  following  well-known  results  [recall  Eg.  (26)]: 


lL  h  (koR»  ylm  <r) 
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In  Eqs.  (36)  and  (37),  the  quantity  C^m  Aim.  is 


'im,  X'm' 
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where  C(L,  i,  it';  M,  m)  are  Clebsch-Gordan  coefficients. 
Eqs.  (38)  and  (39),  we  have  assumed 

PS  *  pa  ' 

which  is  the  case  of  interest  and  kL(x)  are  again  certain 
modified  spherical  Bessel  functions  defined  as 


,  .L+2  (1) 

K  (x )  =  JL  i  h  (xx)  . 

L  2  L 

Making  use  of  the  above  relationships  in  the  expression  for 
Ha  ( R )  ,  we  find 
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If  Eqs.  (43)  and  (10)  are  now  substituted  into  Eq.  (21), 


obtains  a  system  of  linear  homogeneous  equations  for  the 
determination  of  the  coefficients  '•>  .  Namely,  we  have 
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where  the  matrix  Q  is 
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The  quantities  VA  and  are 
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(46) 


(47) 


Here,  R1^  and  (R^)  '  are  the  wave  function  and  its  derivative, 
respectively,  evaluated  at  p  =  Pa .  In  order  to  express  Q  in 
the  form  of  Eq.  (45) ,  it  is  necessary  to  use  the  Wronskian  of 
the  (modified)  spherical  Bessel  functions, 

w  j  (x)  /  (x) 

W  I  (x)  ,  K  (x)  =  -  1  (48) 

l  *■  2 


In  order  for  the  system  of  equations,  Eq.  (44),  to  have  a 
solution , 


This  secular  equation  determines  E  as  a  function  of  k  and 
thereby  defines  the  band  structure. 

For  each  E  that  satisfies  Eq.  (49) ,  one  can  determine  the 
solutions  of  Eq.  (44)  up  to  a  common  unknown  factor. 

This  latter  is  determined  by  the  requirement  that  the  wave 
function  Eq.  (19)  be  normalized.  After  normalization, 

Eq.  (10)  specifies  the  interior  wave  function  and  Eq.  (12) , 
the  exterior.  With  the  band  structure  and  crystalline  wave 
function  specified,  the  total  crystal  energy  can  be  found. 
With  our  present  initial  model,  normalization  of  the  wave 
function  is  not  necessary. 

Experience  indicates  that  the  determination  of  the  band 
energy  requires  the  solution  of  Eq.  (49)  for  E  at  several 
thousand  values  of  k.  However,  because  of  the  complex  tran¬ 
scendental  character  of  the  determinant,  it  must  be  evaluated 
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fl£ 

U 


at  sufficient  E  values  to  permit  interpolation  for  that  value 

of  E  that  satisfies  Eq.  (49)  for  a  given  £  vector.  Since  the 

ct  8 

determinant  involves  all  the  DLM's,  which  depend  both  on  E  and 
k,  it  is  clear  that  the  d£m's  have  to  be  evaluated  at 

many  tens  of  thousands  of  k  values.  As  a  consequence  the 

ct  3 

time  to  calculate  the  infinite  series  can  dominate  the 
time,  and  therefore  the  cost,  of  the  overall  calculation. 

Thus  minimizing  this  time  is  vital. 


Up  to  now,  the  literature  has  emphasized  only  one  approach  to 

ct  3 

improving  the  rate  of  convergence  of  the  DLM's,  namely  the 
application  of  the  Ewald  method.  While  substantial  improve¬ 
ment  is  possible  this  way,  it  suffers  from  the  requirement 
that  for  a  given  accuracy  either  time-consuming  incomplete 
gamma  functions  must  be  evaluated  or  a  large  number  of  terms 
in  the  infinite  series  must  be  calculated. 


In  the  present  work,  what  appears  to  be  an  entirely  new 
approach  is  being  explored.  The  basis  of  this  approach  rests 
on  the  observation  that  d£m  as  given  in  Eq.  (35)  is  independ¬ 
ent  of  R  subject  only  to  the  constraint  Eq.  (29) . 


Thus  we  can  multiply  Eq.  (35)  by  an  arbitrary  function  f(R) 

and  integrate  on  R  up  to  some  upper  limit  constrained  by 

a  8 

Eq.  (29).  Thus,  DLM  can  be  rewritten  as 


4£l  y  ei  (k+G)*  (a8-aa)  YLM  (k+G) 

s  ht  4  -*•  -*•  2  FAC  (L,M) 
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where  we  have  defined  the  following  functions: 


FL(X)  =  J'1^  fL(y)  iL(xy)'  (51) 

0 

[p.  (k  8)1  E  >  0, 

f l (y )  iL(kQRy)  ,  E  <  0,  (52) 

dy  f  (y)  nQ  (kQRy)  ,  E  >  0, 
dy  f  (y)  Kq (kQRy) ,  E  <  0.  (53) 


of  k,  Eq.  (50)  must  be  evaluated  for 
many  values  of  E,  calculation  time  would  be  prohibitive  even 
with  improved  convergence  of  the  sum.  This  problem  can  be 
minimized  by  use  of  the  identity 

N-l 

i  =  i  /  e  \m  +  /  e 

|  k+G  1  -E  [  k+G  |  m=0  \|k+G|  /  \|k+G| 

which  allows  E  to  be  factored  from  the  sum  except  for  the 

last  terra  in  Eq.  (54).  However,  even  for  modest  values  of 

N,  the  last  term  can  be  made  to  converge  extremely  fast.  The 

actual  implementation  of  Eq.  (54)  will  be  discussed  below  in 

ag 

connection  with  the  computer  program  for  calculating  DLM. 

We  now  turn  to  the  choice  of  f(y).  We  require  first  that  f (y) 
be  such  that  the  integrals  over  y  be  analytically  evaluable. 


.  N 


(54) 


k+G I  -E 


o 


l  0 


Since,  for  a  given  value 


If  numerical  integration  were  required,  the  time  for  evalua¬ 
tion  would  be  prohibitive.  A  second  desirable  property  of 
f(y)  is  that  it  vanish  at  y  *  0  and  y  =  1.  This  property 
implies  that 


J  dy  f(y)  j  ( [  k+G  !  Ry)  ~  - — 


for  large  |  k  +  Gl ,  where  n  depends  upon  the  order  of  the  zero 
at  y  =  0  and  y  =  1. 

Although  a  large  number  of  functions  satisfy  these  conditions 
we  have  concentrated  on  the  following  two: 


f  (y)  =  -L-  yL+2  (1-y2 )P, 
1  2V  ul 


f  (y)  =  f-  yL+1  d-y2)L. 

2  Ll 


For  these  functions,  we  find  for  FL(x),  Eq.  (51), 


f  =  f 


F  (x)  .  . , 


f  =  f 


=  =[3Lff)]2/(f)L 


The  results  for  HL  and  H  will  be  discussed  below.  It  is 
clear  from  Eqs.  (58)  and  (59)  that  R  should  be  chosen  as 


large  as  possible  consistent  with  Eq.  (29).  However,  which 
form  provides  the  greater  convergence  improvement  is  more 
subtle.  To  clarify,  one  observes  that  for  small  argument 


j  (x)  ~  2T.  xn  (60) 

n  r ( 2n+2) 

and  for  large  arguments 

j  (x)  ~  i  sin 
n  x 

\ 

Eq.  (60)  points  up  the  known  result  that  the  larger  n  becomes, 
the  larger  x  must  be  before  Eq.  (61)  starts  to  apply.  Thus, 
in  Eq.  (58)  the  larger  n  becomes  the  more  is  the  number  of  J 

G  values  required  before  the  denominator  can  speed  conver-  \ 

gence.  On  the  other  hand,  in  Eq.  (59)  the  order  of  the 
spherical  Bessel  function  is  small  but  the  convergence  factor 
in  the  denominator  is  constrained. 

Taking  account  of  these  conflicting  advantages  suggests  that 
Eq.  (58)  is  probably  to  be  preferred  for  low  values  of  L  and 
modest  values  of  p,  while  Eq.  (59)  is  to  be  preferred  for 
larger  values  of  L. 
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IV.  COMPUTER  PROGRAM  DEVELOPMENT 

In  developing  the  computer  program  a  number  of  assumptions 
have  been  made: 

•  the  maximum  l  of  interest  for  the  wave  function  is 
A  *  3  and,  therefore,  the  maximum  L  of  interest 
for  the  structure  constants  is  L  =  6, 

•  the  maximum  number  of  ions  per  unit  cell  is  2, 
though  this  restriction  may  be  readily  removed, 

•  outward  recurrence  relations  are  acceptable  to 
calculate  jA  ( x )  and  (x)  , 

•  f^(y)  with  n  =  5  is  used  for  L  =  0  -  4  and  f2<y) 
is  used  for  L  *  5  -  6, 

•  the  parameter  N  in  Eq.  (54)  is  taken  to  be  10. 


The  last  three  of  these  assumptions  will  be  discussed  in 
Section  VI. 

Once  the  lattice  has  been  specified  (see  Appendix) ,  approxi¬ 
mately  4000  G  values  are  computed  and  sorted  according  to 
length.  A  shell  structure  in  G  space  is  defined  to  test  for 
convergence  of  the  expansion  coefficients.  At  present,  this 
shell  structure  is  externally  prescribed. 

The  first  step  is  to  calculate  expansion  coefficients  based 
on  Eq.  (54).  We  define  the  two  quantities  BR  and  BI ,  for 
M  >  0 ,  by 
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BR (L,M,I) 

= 

P(L,M)  p  /p  1  qxp"  j  \ 

COS  M  tj> 

BI (L,M, I ) 

«► 

G 

”|  G+k  |  2]  L 

sin  (-M)<}> 

(62) 


where  P(L,M)  are  the  usual  associated  Legendre  polynomials 
and  ij>  is  the  azimuthal  angle  of  the  vector  G+k.  If  there  are 
two  ions  per  unit  cell/  we  also  define  the  four  quantities 
Bl,  B2,  B3  and  B4/  for  M  >  0,  as 


' 

Bl (L,M,I) 

cos  M<(>  cos  X 

B2 (L,M, I ) 

.  =  P(L,M)  p  ,B  iqj.y  1  \ 

sin  (-M)4>  cos  x 

B3(L,M,I) 

*  ri?4.v|2ll  L 

cos  M<J>  sin  x 

B4 (L,M, I ) 

d 

G  |J  G+k  |  J 

sin  ( — M <j) )  sin  x 

4 

where 

x  =  (*2"^)  *  (G+ic)  .  (64) 

As  each  shell  is  calculated,  the  fractional  change  in  the 
quantity 


(BR)2  +  (BI)2 

is  determined  for  those  values  of  L,  M  and  I  considered  in 
the  shell.  If  for  a  given  L,  M  and  I,  the  change  is  less 
than  some  prescribed  amount,  that  L,  M  and  I  is  considered 
to  have  converged  and  will  not  be  calculated  for  the  next 
succeeding  shell.  Consequently,  as  we  go  through  the  shells, 
contributions  are  calculated  for  fewer  and  fewer  expansion 
coefficients  until  all  have  converged.  At  present,  no 
testing  is  done  on  Bl-4. 

The  program  for  this  part  of  the  calculation  has  been  opti¬ 
mized  to  a  large  extent  but  some  improvement  cannot  be  ruled 
out.  Most  of  the  testing  of  this  portion  of  the  program 
assumed  a  fee  lattice  with  a  volume  of  2ir* ,  one  ion  per  unit 
cell  and  a  k  of 
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k  =  (0.3, 0.4, 0.5). 


(65) 
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The  maximum  value  for  the  parameter  R  is  then 

R  =  4.4.  (i 

Q 

The  reason  for  these  choices  is  that  published  values  for 
the  structure  constants  at  the  energy 


E  =  0.656  (67) 

were  available  for  comparison.  We  found  that  with  these 
choices  (particularly  R) ,  all  the  expansion  coefficients  con¬ 
verged  to  an  accuracy  of  10~3  with  about  1000  G  points.  To 
improve  the  accuracy  to  10"4  required  about  2000  G  points. 

If  R  is  decreased  to  2.7,  some  of  the  expansion  coefficients 
did  not  converge  to  an  accuracy  of  10”3  even  with  4000  G 
points.  This  point  will  be  discussed  in  Section  V. 


m 


Once  the  expansion  coefficients  for  a  given  k  value  have  been 
evaluated,  the  determinant  for  a  predetermined  set  of  energies 
can  be  calculated.  To  do  so  a  set  of  coefficients  CR(L,  M) , 

Cl (L,  M)  and  if  appropriate,  C1-4(L,  M)  are  calculated.  These 
are  defined,  as  in  Eqs.  (62)  and  (63),  respectively,  by  the 
substitution 


[l^l2]1  *  [l  S+K  [ 2] 


The  calculation  is  done  for  the  first  shell  only  and  no  test 
of  convergence  is  necessary  because  of  the  rapidity  of  conver¬ 
gence.  Based  on  Eq.  (54),  the  appropriate  combinations  are 


u 


DX (L,M) 


_  4 ff  y 

Q 


1-1 

E  BX(L,M,I) 


CX (L,M) 


where  X  =  R,  I ,  1-4.  Simple  linear  combinations  of  these 
quantities,  divided  by  HL,  then  yield  DLM,  Eq.  (50),  except 
for  the  H  term.  We  note  that  for  E  <  0,  HL  is  simply  (x  =  kQR) 


f  =  f,  :  H  =  I 


L  L+u  +  1 


P  +  l 
(x)/x  , 


t  =  £*  :  Hr.  =  [rr.(f)]2/(f)L- 


For  L  =  0,  with  fjjy)  given  by  Eq.  (56)  and  ii  =  5,  the  calcu¬ 
lation  of  H,  Eq.  (53),  is  trivial. 


The  portion  of  the  program  for  calculating  the  structure 
constants  has  been  tested  and  improved  using  the  standard 
lattice  defined  above.  Typically,  the  structure  constants 
have  the  same  relative  accuracy  as  the  expansion  coefficients, 
The  only  exceptions  may  be  structure  constants  with  abnor¬ 
mally  small  values.  In  Table  1,  we  present  the  structure 
constants  as  given  in  Ref.  8  and  our  calculation  of  them 
using  error  measures  of  10"4  and  10"3,  respectively.  Since 
the  normalization  in  Ref.  8  is  different  from  ours,  Eq.  (50), 
we  have  converted  our  normalization  for  ease  of  comparison. 
Notice  that  for  M  >  0  (M  <  0) ,  the  structure  function  is 
related  to  BR  (BI),  Eq.  (62). 

Once  the  structure  constants  have  been  determined,  the  matrix 
elements  readily  follow.  The  quantities 
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TABLE  1.  STRUCTURE  CONSTANTS,  fee,  VOLUME  »  2H5 
it  ■  (0.3,  0.4,  0.5),  E  *  0.656 


Djj^tOavl  s ) 

Dm(i°"4) 

°lm( 1°“3> 

0.3637069 

0.363707 

0.363707 

0.4729255 

0.472924 

0.472922 

0.5585931 

0.558592 

0.558583 

0.3855073 

0.385506 

0.385504 

0.1383943 

0.138390 

0.138385 

0.3854948 

0.385489 

0.385477 

0.2294027 

0.229400 

0.229400 

0.2425004 

0.242496 

0.242488 

-0.1231771 

-0.123176 

-0.123175 

0.3342447 

0.334248 

0.334245 

0.7787625 

0.778762 

0.778763 

0.3760202 

0.376019 

0.376012 

-0.3433054 

-0.343305 

-0.343308 

0.3140171 

0.314016 

0.314011 

-0.0644856 

-0.064484 

-0.064482 

-0.3585324 

-0.358534 

-0.358529 

-0.1350521 

-0.135051 

-0.135050 

-0.5645986 

-0.564611 

-0.564689 

0.0665038 

0.066486 

0.066451 

0.4518756 

0.451877 

0.451905 

-0.2631646 

-0.263159 

-0.263139 

0.4007033 

0.400704 

0.400723 

-0.1747519 

-0.174749 

-0.174753 

0.1291246 

0.129140 

0.129196 

0.2032239 

0.203229 

0.203258 

0.49  58597 

0.495869 

0.495955 

-0.0768520 

-0.076852 

-0.076854 

0.8811268 

0.881139 

0.881165 

w 


STRUCTURE  CONSTANTS,  fee,  VOLUME  -  2n3 
k  »  (0.3,  0.4,  0.5),  E  »  0.656  (CONCLUOEO) 


Dlh(DovI  s  ) 

K> 

1 

O 

5 

O 

0.2049815 

0.204982 

0.204986 

-0.9030668 

-0.903089 

-0.903093 

0.6678767 

0.667894 

0.668039 

-0.6314083 

-0.631422 

-0.631417 

-0.4744635 

-0.474484 

-0.474487 

-0.8026519 

-0.802659 

-0.802677 

-1.5614385 

-1.56146 

-1.56226 

0.0352270 

0.035233 

0.035293 

4.3426328 

4.34269 

4.34528 

-2.2561097 

-2.25614 

-2.25619 

-0.3471268 

-0.347127 

-0.347125 

4.5963681 

4.59643 

4.59700 

-2.3159161 

-2.31596 

-2.31813 

3.5271685 

3.52722 

3.52776 

-1.1349680 

-1.13498 

-1.13508 

3.0624581 

3.06251 

3.06295 

0.1585912 

0.158593 

0. 1 58807 

-4.1491049 

-4.14916 

-4.14964 

0.9372909 

0.937309 

0.937489 

-1.4263739 

-1.42640 

-1.42644 

0.3302557 

0.330258 

0.330567 

CB(2m,A'm'  ,L)  =  FAC  (L ,  M) 


(71) 


are  calculated  early  in  the  program  and  stored  for  use  as 
necessary.  For  a  prescribed  set  of  energies,  the  quantities 
v£,  Eq.  (46),  and  U®,  Eq.  (47),  are  also  calculated  and 
stored.  The  matrix  elements  Q,  Eq.  (45),  are  computed  and 
the  determinant  calculated.  This  final  process  has  been  pro¬ 
grammed  but  not  checked  extensively. 

Finally,  a  word  about  V“  and  u£  is  required.  An  existing 
Thomas-Fermi  program  provides  the  basic  information  ior  the 
wave  function  and  its  derivative,  namely  E  (TF ) ,  PSI(fc,a)  and 
DPSI/DR ( i,  a)  which  are  related  to  the  above  quantities  by 


$ 


I 


E  =  -E (TF)  , 

Rj  *  PSI  U,a)  , 

(R?)  '  =  DPSI/DR  (*, a)  -  -L  PSIU,a).  (72) 

P<x 


3 

£ 


51 


& 


The  calculation  of  v£  and  U®  using  the  provided  quantities 
has  been  programmed  but  not  checked. 

The  complete  program  was  run  using  carbon  wave  functions  with 
p0  *  1.26  and  a  set  of  energies  from  -4  to  +2  at  intervals 
of  0.1.  A  simple  cubic  lattice  was  assumed  with  a  volume  of 
19.2  and  the  value  of  R  was  taken  to  be  2.5,  near  its  maximum 
value  for  this  lattice.  The  k  vector  was  chosen  to  be 


k  =  (0.1, 0.2, 0.3)  .  (73) 

The  determinant  as  a  function  of  energy  implied  a  number  of 
zeros  in  this  energy  range  but  the  actual  zeroes  were  not 
evaluated . 
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IMPROVED  STRUCTURE  CONSTANT  EVALUATION 


The  numerical  work  carried  out  so  far  strongly  suggests  that 
the  approach  using  Eq.  (50)  is  adequate  for  one  ion  per  unit 
cell.  However,  for  crystal  structures  with  more  than  one  ion 
per  unit  cell  the  maximum  value  of  R  is  in  general  decreased. 
For  example,  if  in  our  test  fee  lattice,  we  considered  the 
diamond  structure,  the  maximum  value  of  R  is  2.7.  Calcula¬ 
tions  with  this  value  of  R  indicated  poor  convergence  prop¬ 
erties  for  the  expansion  coefficients.  The  problem  arises, 
of  course,  for  small  values  of  I  in  Eqs.  (62)  and  (63) . 

A  possible  solution  to  this  problem  as  well  as  an  improved 
one  ion  per  unit  cell  calculation  is  based  on  a  slightly 
altered  expansion  [cf.  Eq.  (54)], 


N-l 

— i —  =  y  (E+A2)m _ i _ 

I  k+G  I  2  -E  m=0  [|  k+G  I  2+X2]m+1 

(E+X2)N _ 1 _ f  (74) 

I  k+G  I  2-E  [l  k+G  I  2+X2]  N 

where  X  is  arbitrary.  We  define  the  expansion  coefficients 
BR,  BI  and  Bl-4  just  as  in  Eqs.  (62)  and  (63),  respectively, 
with*  the  substitution 


[< 


(75) 


Corresponding  changes  also  occur  in  the  CX  coefficients, 
Eq.  (68)  . 


For  the  larger  values  of  m,  we  use  exactly  the  same  techniques 
as  discussed  above  since  convergence  should  be  reasonably 


rapid  as  long  as  X  is  not  large  compared  to  tha  lowest  values 
2 

of  G  .  For  the  first  few  m  values,  say  m  =  0,  1,  2,  we  pro¬ 
pose  to  calculate  the  expansion  coefficient  in  coordinate 
space.  We  observe  first  that 


[i  s+a  i2+x2]: 


1  _d_  1 

2X  dx  |£+g|2+x2 


1 

-  1  1 

d2 

_  i 

d 

1 

[ik+G  1  2+X2  ]  3 

ICM 

*< 

I  00 

CM 

_ i 

X 

dX 

Ik+G 1 Z+XZ 

Next  we  equate  the  two  representations  for  H ( R) ,  Eqs.  (33) 

2 

and  (22)  [for  negative  energy  E  =  -X  ] 

4  Jr  V  i  (k+G)  *R 


4jr  y  e1(ittt3)  iara 
♦  |  k+G I z  +  Xz 


>  *  >  > 
— X I RM+a„— aQ-R 


_  y  eik*RN  e~A  I  RN't'aa~aB 


WV*1 


and  make  a  YLM(R)  decomposition,  leading  to 


r"  i  (k+G)  •  (aft-a  ) 

2 ^  - - - - - — —  jL  ( I  k+G  I  R)  Y*M(k+G) 

♦  I  k+G I  2  +  X2 


a  e 

yii  iR 

’  5a6  6L0 

— ,  f 

♦  f 

“-L  I 

eik,RN  ' 

N 
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where  the  prime  on  the  sum  means  the  N  =  0  term  is  deleted 
if  a  =  6  and  the  vector  V  is 


V  =  R^j  +  aa  -  a^-  (8  0) 

Multiplying  by  f(R)  and  integrating  on  R,  we  find  an  alter¬ 
nate  expression  for  BR (L ,  M,  1),  namely. 


BR (L,M, 1)  =  JL 
Air 


X  H  6ag  6L0  SM0 


(81) 


+  X  ^  P(L,M)  Hl  Kl(XV)  Re 
N 

-V  t  > 

where  <}>  is  the  azimuthal  angle  of  the  vector  V,  HL  is  given 
in  Eq.  (70)  with  x  =  XV  and  H  is  the  integral  given  in 
Eq.  (53)  for  E  <  0  with  kQ  *  X.  The  expression  for  BI(L,M,1) 
is  obvious,  based  on  Eq.  (81).  Bl-4(L,M,1)  are  somewhat  more 

complicated  since  they  require  the  calculation  in  coordinate 
space  of  both  a  *  1,  8=2  and  a  =  2,  B  =  1.  Once  the  expres¬ 
sions  for  1=1  are  derived,  the  results  for  1=2  and  1=3 
follow  from  Eqs.  (76)  and  (77),  respectively.  These  calcula- 

<v 

lations  are  straightforward  since  the  derivatives  of  IL  and 
can  be  expressed  in  terms  of  modified  spherical  Bessel 
functions . 


e1K  KN  e 


The  advantage  of  this  scheme  results  from  the  fact  that  KL(x) 
falls  off  exponentially  for  large  x.  A  few  calculations  for 
selected  values  of  L  and  M  have  been  carried  out  with  very 
encouraging  results.  However,  no  detailed  examination  of  this 
scheme  has  yet  been  undertaken. 
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CONCLUSIONS 


What  remains  to  be  done  to  complete  the  initial  model  is 
(1)  to  complete  and  thoroughly  check  out  the  computer  program 
for  computing  E  vs  k;  (2)  to  improve  the  calculation  of  jL(x) 
and  lL(x)  *or  small  arguments;  (3)  to  investigate  further  the 
choice  for  y  [Eq.  (56)]  and  the  split  between  f]_  and  f2  as 
a  function  of  L;  (4)  to  determine  the  sensitivity  of  E  vs  k 
as  a  function  of  errors  in  the  structure  constants  and  the 
wave  functions;  (5)  to  investigate  extensively  the  scheme 
discussed  in  Section  V;  and  (6)  to  sum  over  the  bands  to 
calculate  the  total  band  energy. 

Once  the  initial  model  is  completed,  tests  of  its  usefulness 
are  needed,  namely:  (1)  the  calculation  for  a  given  volume 
and  lattice  type  must  have  a  reasonable  cost;  (2)  for  a  vari¬ 
ety  of  known  substances,  the  model  must  correctly  predict  the 
stable  crystal  structure  (note  that  this  does  not  mean  the 
total  cohesive  energy  need  be  calculated  accurately,  only  that 
the  observed  structure  has  the  lowest  energy  according  to  our 
model);  (3)  for  known  phase  transitions  at  elevated  pressures, 
we  should  identify  the  fact  that  a  phase  transition  occurs. 
Note  that  this  does  not  mean  we  could  accurately  predict  at 
what  pressure  the  transition  occurs — but  only  estimate  it 
approximately. 

In  applying  these  tests,  we  should  be  willing  to  be  wrong 
in  some  fraction  of  the  cases,  since  the  model  is  to  serve 
only  as  a  guide  for  experiments.  However,  if  the  model 
failed  most  of  the  tests,  it  will  probably  be  necessary  to 
consider  some  or  all  of  the  ignored  energy  contributions, 
and  possibly  even  an  improved  ion  potential. 


APPENDIX 


A  crystal  lattice  is  defined  by  its  three  primitive  basis 
vectors,  a,  S  and  c.  The  possible  relationships  between 
these  vectors  fall  into  14  categories,  called  Bravais  lat¬ 
tices.  As  input  to  the  computer  program  one  prescribes  the 
lattice  type,  the  volume  of  a  unit  cell,  ft,  and  any  needed 
length  ratios.  From  this  information,  the  primitive  lattice 
vectors  are  calculated  and  from  these  the  reciprocal  lattice 
vectors,  a*,  b*  and  c* , 


=  If  b  x  ?,  b*  =  ll  ?  x  a,  ?*  =  If  a  x  b. 


A  momentum  space  vector,  G,  is  then  given  by 


G  *  n^  a*  +  n2  b*  +  n3  c*. 


while  a  direct  lattice  vector,  RN ,  is  given  by 


=  Nj^a  +  N2  b  +  N3  c. 


Here,  j  n  ,  n2,  n3J  and  j  ,  n2,  N3  j  are  sets  of  three  inte¬ 
gers.  For  each  ion  in  the  unit  cell,  its  charge  (za) , 
radius  (pa)  and  position  (aa)  are  specified.  Here,  a  *  1, 
...N,  with  N  being  the  total  number  (<2)  of  ions  in  the  unit 
cell . 
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